Imposing Field Boundary Conditions in Md Simulation of Fluids: Optimal Particle Controller and Buffer Zone Feedback

We formulate a method for imposing continuous field boundary conditions on an MD simulation with little or no disturbance to its dynamics. Our approach combines the previously developed thermodynamic field estimator, which extracts macroscopic fields from particle data, with a novel procedure, optimal in the sense of least disturbance, for imposing prescribed continuum field boundary conditions on the atomistic system. By means of feedback control and assuming molecular chaos in fluids, that disturbance can further be eliminated entirely, thus providing an exact solution for general steady-state fluid problems, where the desired fields at the real boundary is achieved by adjusting actions in a region separated from the real boundary by a buffer zone. OPTIMAL PARTICLE CONTROLLER (OPC) A basic problem in fluid simulation, which was first investigated by O’Connell and Thompson[1], then Hadjiconstantinou and Patera[3], and then by us[4, 5], is the following: we would like the particles in an MD simulation to satisfy the following distribution in their positions and velocities fxi;vig, dP = f x;vj (x); T (x); v̄(x) dxdv = (x)dx (2 T (x))3=2 exp jv v̄(x)j2 2T (x) ! dv+f (2)dxdv; (1) at a prescribed boundary @C, where (x), T (x), v̄(x) are the macroscopic density, velocity and temperature fields. For this to be sensible, @C should be in small-gradient regions where (x) etc. are well-defined; conversely, if that is true, the transition from particle to continuum description like the Navier-Stokes equation, can be shown to be exact by the Chapman-Enskog expansion, where f (2) is the the second-order correction to the local Maxwellian that contains the spatial gradients of these fields. So, both continuum and particle descriptions should be valid when close to @C, and (1) is the link that allows us to couple the two. In a previous paper[4] we have shown how to solve the opposite problem, that is how to infer the macroscopic fields in C from current particle data, using a technique named thermodynamic field estimator (TFE); the upshot is that we can accurately determine the macroscopic fields at an inner boundary @U inside C in real time even when the particle data contains considerable “thermal noise”. And so @U provides the BC for a continuum solver in the classical domain decomposition and alternating Schwarz iteration formalism[2], where @C @U is the continuum-MD overlap region. Because of space limitation we will not further discuss the TFE here, nor the full Schwarz iteration implementation which we have carried out for Couette flow, but will focus on just the second part of the iteration, which is how to impose continuous field BC on C. Interested readers should follow the pioneering paper of Hadjiconstantinou and Patera[3], and much detailed reports of our work in [4, 5]. Suffice it here to note that in order for the whole scheme to work to second order accuracy in (1), we can simply use the first term in (1) for TFE, ignoring f (2). Now, suppose one has inferred the current fields to be 0(x), T 0(x), v̄0(x) on @C using TFE, but wants the fields to be (x), T (x), v̄(x), how should one modify the particle coordinates such that the desired distribution is achieved? That is, say there is a random variable set fXng satisfying distribution f(X), but we want them to satisfy distribution g(X), so we begin to replace Xn’s by Yn’s hoping that fYng will satisfy g(Y ), what should be the optimal Xn ! Yn transformation? One may ignore Xn and use only g(Y ) to resample the random variable Yn. While this would indeed lead to a BC imposing scheme, it is quite conceivable that ignoring the current state will result in a procedure that strongly disturbs the particle dynamics. If decoupling the two distributions, or equivalently, the two sets of particle coordinates, is not a good idea, then one should look for a way to relate them such that the disturbance is minimized. We now introduce a quantitative measure of this disturbance: B = Xn j∆vnj2 =Xn jvout n vin n j2; (2) where ∆vn is the change in particle velocity of the nth particle (vin n and vout n are the velocity before and after the transformation). This quantity characterized the artificial disturbance to particle dynamics which should be minimized as much as possible. As an example let us consider a simple problem of 1D heat conduction, for which we wish to impose a boundary condition of high temperature Th at x = 0 and low temperature Tl at x = 1. Intuitively one could imagine doing the following: when a particle crosses either boundary x = 0 or 1, give it a random velocity drawn from distribution (1) with parameters v̄ = 0 and T = Th or Tl. However, this procedure does not work. When implemented in practice in the case of homogeneous heating (Th = Tl), results show that the bulk temperature reaches a value of Th=2. The scheme fails because one is dealing with conditional probability. The speed distribution of atoms which cross the boundary is different from the speed distribution of atoms in the bulk. Instead it is weighted by the normal velocity, dP = v exp v2 2T dv R+1 0 v exp v2 2T dv = v T exp v2 2T ! dv; 0 < v < +1: (3) hv2i from distribution (3) is 2T , not T . Thus, if we sample the boundary crossing atoms using bulk distribution (1) with parameter T = Th = Tl, the energy can only be balanced in a statistical sense when the bulk temperature reaches Th=2. A more subtle defect of this scheme is that, for whatever the incoming velocity vin of the particle before hitting the boundary, a new velocity vout is drawn from a given distribution, say g(v), entirely independent of vin. Thus if we evaluate the disturbance to particle dynamics using (2), it is always substantial. In the case of homogeneous heating, even when the correct distribution (3) is used in drawing vout’s and the system has reached the desired temperature T = Th = Tl, the scheme continues to disturb the particles by giving each boundary crossing atom a new vout. On the other hand, if we just let vout = vin, i.e., do nothing, the system temperature stays at Th! A more intelligent particle controller should automatically tune down its influence as the system approaches the desired state, a behavior we may call the coalescence property. We now formulate the above ideas mathematically. Suppose fXng conform to distribution function f(X): dP ( < X < + d ) = f( )d ; (4) and we would like the series to conform to a different distribution g. We propose to achieve this by replacing every X with another number Y which is distributed according to g. Thus the goal is to find a transformation T , T : Xn ! Yn; (5) with the requirement that if fXng conforms to distribution f(X), fYng will conform to g(Y ): dP ( < Y < + d ) = g( )d : (6) There are many possible T ’s. However our previous discussion shows that the following property is desirable: if f g, then T gives Yn Xn. Generally, to incorporate the idea of minimally disturbing the dynamics, we adopt the reasonable criterion that B[T ] = h(Y X)2i (7) be minimized among all possible T ’s. We see that if T is randomly drawing Y from g(Y ) without referencing to X , it satisfies the basic requirement (6) but not the coalescence property. We call this T1 transformation. To incorporate coalescence one may consider T2, T2 : 8<: Y = X : p g(X) Kf(X) draw Y randomly from g(Y ) : p > g(X) Kf(X) where p is a random number uniformly distributed over [0; 1], and K is a constant (scheme fails if K does not exist) such that Kf(X) g(X) for 1 < X < +1: (8) But is T2 the best one? We propose the following transformation, in the form of an implicit relation, T3(X ! Y ) : Z X 1 f( )d = Z Y 1 g( )d : (9) It can be checked that Yn indeed conforms to distribution g if the incoming random variable Xn conforms to f , and it satisfies the coalescence property. We believe that T3 is the mathematically optimal transformation which minimizes (7). A tentative proof is given in [5]. This assertion should come as no surprise because it is the only one-to-one continuous mapping which satisfies (6) without extra randomness like in T2. Our experience shows that while T3 works rather well, it is an implicit algorithm and could be computationally demanding. One may choose to operate T in a finite-volume region or just at a certain boundary, which are called bulk or boundary particle controllers respectively. The subtleties of boundary particle controllers are already mentioned in the first example. When the average disturbance to microscopic particle dynamics, defined by some reasonable measure, due to the operation of a family of particle controllers serving certain fixed purpose, is minimized, it is called the optimal particle controller (OPC). We believe that for oneor decoupled multi-dimensional distributions, T3 is OPC. The explicit formulas for general 3D T (x), v̄(x) boundary particle controller are given in [5], while (x) coupling involves particle insertion or removal which is more difficult to treat (see [4]). But explicit (x) boundary condition is rare. Also, those derivations ignore f (2), which cause certain problems, and we will derive more accurate versions by factoring f (2) into the distribution function. In all cases the general discussions about OPC hold. EXTENDED BOUNDARY CONDITION In last Section a method was proposed to control the field boundary condition of an MD simulation that results in least disturbance to the particle dynamics in the sense of (2). Nevertheless the disturbance still exists for particles in the skin region near the boundary, which can be shown to be proportional to the rate of dissipation in the system. In this Section we will formulate a procedure which eliminates that disturbance entirely at the specified boundary, thus providing an exact theoretical solution to the problem. Since any artificial action necessarily alters the particle dynamics in the vicinity where it is imposed, the best one can do is to act some distance away from the intended boundary and cause the macroscopic field at the boundary to be what is prescribed. This can be done using a three-region approach which we will call the Extended Boundary Condition (EBC), shown in Fig.1, and through a feedback control mechanism. The physical region of interest, C (core), is surrounded by a buffer zone B. Artificial actions are applied on an outer MD region A which is sufficiently separated from C, its aim being to induce the prescribed field boundary conditions on the core boundary @C. Due to molecular chaos in fluids, short wavelength Control Algorithm f g f Field Estimator Particle Controller A ct io n R eg io n( A ) B uf fe r R eg io n( B )